Exploring the most promising anti ‐ Depressant drug targeting Microtubule Affinity Receptor Kinase 4 involved in Alzheimer’s Disease through molecular docking and molecular dynamics simulation

Alzheimer’s Disease (AD) is the prevailing type of neurodegenerative illness, characterised by the accumulation of amyloid beta plaques. The symptoms associated with AD are memory loss, emotional variability, and a decline in cognitive functioning. To date, the pharmaceuticals currently accessible in the marketplace are limited to symptom management. According to several research, antidepressants have demonstrated potential efficacy in the management of AD. In this particular investigation, a total of 24 anti-depressant medications were selected as ligands, while the Microtubule Affinity Receptor Kinase 4 (MARK4) protein was chosen as the focal point of our study. The selection of MARK4 was based on its known involvement in the advancement of AD and other types of malignancies, rendering it a highly prospective target for therapeutic interventions. The initial step involved doing ADMET analysis, which was subsequently followed by molecular docking of 24 drugs. This was succeeded by molecular dynamics simulation and molecular mechanics generalised Born surface area (MMGBSA) calculations. Upon conducting molecular docking experiments, it has been determined that the binding affinities observed fall within the range of -5.5 kcal/mol to -9.0 kcal/mol. In this study, we selected six anti-depressant compounds (CID ID ‐ 4184, 2771, 4205, 5533, 4543, and 2160) based on their binding affinities, which were determined to be -9.0, -8.7, -8.4, -8.3, -8.2, and -8.2, respectively. Molecular dynamics simulations were conducted for all six drugs, with donepezil serving as the control drug. Various analyses were performed, including basic analysis and post-trajectory analysis such as free energy landscape (FEL), polarizable continuum model (PCM), and MMGBSA calculations. Based on the findings from molecular dynamics simulations and the MMGBSA analysis, it can be inferred that citalopram and mirtazapine exhibit considerable potential as anti-depressant agents. Consequently, these compounds warrant further investigation through in vitro and in vivo investigations in the context of treating AD.


Introduction
Alzheimer's disease (AD) is a rapidly expanding neurodegenerative condition characterised by progressive and irreversible deterioration of the central nervous system, resulting in cognitive and behavioural deficits [1][2][3].The process of diagnosing a condition relies on the analysis of historical case data, clinical evaluations, neuroimaging techniques, and blood analyses.The brain has two primary categories of aberrant proteins, namely amyloid-beta (Aβ) and tau.The disease advances as a consequence of the accumulation of these proteins [4][5][6][7].At now, the therapeutic options for AD are limited to a total of five drugs.Among them, four belong to the class of cholinesterase inhibitors, while the last prescription, memantine, functions by inhibiting glutamate formation.Cholinesterase inhibitors are administered to patients diagnosed with mild to moderate AD, while memantine is prescribed for individuals in the mild to severe stages of AD [8][9][10].
In the present work, the researchers examined the possibility of microtubule affinity regulation kinase 4 (MARK4) as a target for AD.Previous studies have indicated an observed overexpression of MARK 4 in AD, and this overexpression has been linked to the progression of AD.The discovery of the crystal structure of MARK4 was facilitated by the utilisation of the pyrazolopyrimidine inhibitor [11].The overexpression of MARK4 has been found to be accountable for the phosphorylation of the tau protein at the Ser262 location, a crucial step in the binding of microtubules to tau proteins [11].In their study, Lund et al. 2014 observed the co-localization of phosphorylated MARK4 and phosphor-tau Ser262 in granulomatous (GVD) formations found in samples of AD [12].The findings of this study proposed the development of a very effective MARK4 inhibitor for the treatment of AD [12].In the context of this study, it is posited that MARK4 represents the most suitable candidate for addressing neurodegeneration in AD.
AD is widely regarded as a condition that lacks a definitive cure.To date, a total of 24 selective serotonin reuptake inhibitors (SSRIs) have been examined in the existing literature as potential treatments for some symptoms associated with major depressive disorder [13].The development of AD is frequently associated with the early manifestation of depression, characterised by a chronic relapsing-remitting pattern that tends to worsen and become more extended as individuals grow older.This progression significantly heightens the probability of acquiring AD.Late-life depression has been found to be associated with reduced levels of neurotrophins, particularly brain-derived neurotrophic factor.Additionally, it is characterised by the activation of neuroinflammatory pathways, heightened secretion of pro-inflammatory cytokines and C-reactive protein, and an increased cortical amyloid burden.These factors collectively contribute to the process of neurodegeneration.Historically, there was a prevailing belief that the sole determinant of the symptoms associated with AD was cholinergic dysfunction, as evidenced by references [14][15][16].Geeldenhuys, Van der Schyf, and Ramirez have highlighted the significance of the 5-Hydroxytryptamine (5-HT)1, 5-HT4, 5-HT6, and 5-HT7 receptor classes in the context of cognitive enhancement.In comparison to conventional antidepressants such as tricyclics, selective serotonin reuptake inhibitors (SSRIs) exhibit a more desirable side effect profile and have been authorised as primary therapeutic interventions for depressive disorders [14][15][16].There exists compelling data indicating that SSRIs have demonstrated a beneficial impact in delaying the onset of AD in those with a prior history of depression.In addition to the reduction of serotonin transporter function, newly developed serotonergic antidepressants such as vortioxetine exert an impact on serotonin receptors, specifically through the antagonism of 5-HT7 receptors.Vortioxetine significantly improved cognition when compared to other traditional anti-depressants in depressed adults with moderate AD when compared to placebo [17].Based on the recent meta-analysis finding, sertraline and mirtazapine can be considered as an alternative treatment for depression in AD [18].Even the combination of antipsychotics (such as risperidone and quetiapine) and mirtazapine can be used for the management of AD [19,20].
Upon reviewing the previous successful results of SSRIs in managing AD symptoms, the present study conducted molecular docking analyses for a total of 24 antidepressant compounds.Subsequently, utilising the top five binding affinities, molecular dynamics simulations were performed for all six SSRI drugs (CID ID -4184, 2771, 4205, 5533, 4543, and 2160).Both CID 2160 and CID 4543 were subjected to simulation as their binding affinities were found to be same.CID ID-3152 (Donepezil), an authorised medicine for AD, was selected as the control in our study.The primary objective of this study is to identify the most efficacious SSRI medication that specifically targets the MARK4 protein implicated in AD.

Selection of SSRI drugs
We took 24 anti-depressant drugs from the literature based on their mode of action as well as easy availability and well establishment in the global market and it is shown in Table 1.

ADMET analysis of 24 SSRI drugs
The pharmacokinetics and Lipinski's rule are central to accelerating drug discovery and development.As 24 SSRI drugs are already in the market, but still we performed ADMET analysis with the SWISS ADME server which is shown in S1 Table and we even plotted the radar graph for the the binding affinity of the top six drugs and Donepezil.

Molecular docking of 24 SSRI drugs
Molecular docking was performed with AutoDock 4.2.6.First, we converted PDB structures to PDBQT format, added Gasteiger charges to the ligands, and then loaded MARK4-associated protein (5ES1) to AutoDock MGL tools.A grid box with x, y, and z directions was generated to cover the active sites of the target protein, keeping enough space for the rotational and translational movement in the ligands and keeping the grid spacing at 0.375 Å, and flexible docking was performed.The outputs were analysed and visualised in AutoDock 4.2.6.

Molecular dynamics simulation of 6 SSRI drugs
Desmond, a program from Schro ¨dinger LLC [21], was used to assess the effect of the solvation model on the binding affinity between the best-screened ligands and target protein.Molecular docking is a static representation of interactions between protein and ligand in a vacuum [22].In contrast, molecular dynamics simulations use Newton's classical equation of motion to compute atom motions over time.An attempt was made to assess the binding status of all complexes with the target protein (MARK4; PDB ID:5ES1) in the physiological environment predicted using simulations [23,24].The target protein complexes with the drugs were subjected to molecular dynamics simulations for 100 ns each.Pre-processing of the protein-ligand complexes was done using the Protein Preparation Wizard of Maestro software to optimise and minimise the structure of the protein-ligand complexes.All of the systems were prepared using the System Builder tool.TIP3P water (Transferable Intermolecular Interaction Potential 3 Points) was chosen as a solvent model with an orthorhombic box.The simulation was run using the OPLS 2005 force field [25].To simulate physiological conditions, 0.15 M sodium chloride (NaCl) was added to the solvated complexes to neutralise the system.The NPT ensemble was used throughout the simulation, with a temperature of 300 K and a pressure of 1 atm.

Molecular Mechanics Generalised Born Surface Area (MMGBSA) calculations
The binding free energies (Gobind) of the anchored complexes were calculated using the first Molecular Mechanics Generalized Born Surface Area (MMGBSA) module (Schrodinger suite, LLC, New York).Binding free energies were calculated using the OPLS 2005 force field, VSGB solvent model and rotamer search method [26].After performing MD, 10 sec was used to select the MD orbital frame.The total free binding energy is calculated using the below-mentioned formula: ΔGobind = binding free energy, Gcomplex = free energy of the complex, G-protein = free energy of the target protein, and Gligand = free energy of the ligand.
Ethical Statements: This article does not contain any studies involving human or animal participants performed by any of the authors.

Pharmacokinetics and drug likeliness profile of 24 SSRI drugs
Lipinski's is one of the most essential drug discovery and development factors.In this study all the 24 drugs are available in the market, which means drugs follow Lipinski's rule.In this study, we generated radar plots for the binding affinity of the top six drugs (CID ID -4184, 2771, 4205, 5533, 4543, and 2160) lies between -9 to -8.2 kcal/mol as well as radar plot of Donepezil which is already approved for AD.We use the SWISS-ADME web server for which the requirement is the SMILES notation of the selected compounds, which is used as an input based on the database predicting their pharmacokinetic property and plotted the results as a radar plot through which we can analyse the compound's lipophilic property, size, polarity, insolubility, instauration and flex properties.The more it is towards the centre, the more it is considered a good result.The radar plots are shown as Figs 1-7.

Molecular docking analysis of 24 SSRI drugs
First, molecular docking was conducted to validate the intermolecular interactions between the MARK4 protein and 24 drugs.The docking results were arranged in the order of lowest binding energy with each target protein, as shown in Table 2.The binding affinity of the top six drugs lies between -9 to -8.2 kcal/mol.After analysing the binding modes for the top six protein-ligand complexes, 5SE1-2160 formed an h-bond interaction with ASP196 with a bond length of 2.8 Å. 5SE1-2771 formed an h-bond interaction with ILE62 with a bond length of 2.8 Å. 5SE1-3152 formed an h-bond interaction with GLY65 and ILE62 with bond lengths of 4.0    8) 4184 ( 9) 2771 ( 10) 4205 ( 11) 5533 ( 12) 4543 ( 13) 2160 ( 14) 3152. https://doi.org/10.1371/journal.pone.0301179.g008 Å and 3.9 Å, respectively.5SE1-4184 formed a bond length with ASN183, which was 4.3 Å. 5SE1-4205 formed an h-bond interaction with GLU182 ASN183 with bond lengths of 2.8 Å and 3.3 Å, respectively.5SE1-5533 formed an h-bond interaction with ASP196 with a bond length of 2.7 Å. The compound 4553 showed no h-bond interaction with the active binding site of the protein.Besides this h-bond interaction, all the compounds show salt bridges, as shown in the ligand interaction figures.Moreover, the compound 5533 additionally showed halogen interaction.These six drugs are considered further for molecular dynamics simulation, and we considered donepezil as a control, as it is an established drug for AD.The conformation of docked Sertraline, Fluoxetine, Escitalopram, Fluvoxamine, Paroxetine and Citalopram with the MARK4 is shown in the Figs 8-14.The suitable docked pose was selected by mimicking the crystal structure of the MARK4 complex with pyrazolopyrimidine inhibitor  8) 4184 ( 9) 2771 ( 10) 4205 ( 11) 5533 ( 12) 4543 ( 13) 2160 ( 14) 3152.8) 4184 ( 9) 2771 ( 10) 4205 ( 11) 5533 ( 12) 4543 ( 13) 2160 ( 14) 3152. https://doi.org/10.1371/journal.pone.0301179.g011 pyrazolopyrimidine inhibitor is interacting with.The analysis of docked conformations clearly indicates that all six drugs bind deeper into the cavity, and perhaps decrease the accessibility of MARK4 which may be responsible modulation of its biological functions.

Molecular dynamics of 24 SSRI drugs
Based on the top six binding free energy of Sertraline, Fluoxetine, Escitalopram, Fluvoxamine, Paroxetine and Citalopram to the MARK4 which was found to be -9 kcal/mol,-8.7kcal/mol,-8.4kcal/mol, -8.3kcal/mol, -8.2 kcal/mol and -8.2 kcal/mol respectively from 24 drugs, we proceed towards further analysis which include molecular dynamics simulations at 100ns followed by MMGBSA.A molecular dynamics simulation was performed to assess the dynamic behaviour of the protein.The simulation system's stability was assessed by examining the Root  8) 4184 ( 9) 2771 ( 10) 4205 ( 11) 5533 ( 12) 4543 ( 13) 2160 ( 14) 3152.Mean Square Deviations (RMSD) values throughout the trajectory.The RMSD number calculates the difference between the two structures and shows how stable the simulated system is.
A lower RMSD number implies a more minor difference between the two systems.Additionally, a significant factor for examining the dynamics of the protein-ligand complexes is the stability of the simulated system, which is indicated by a less prominent variation in the RMSD value.The kinetics and the binding interactions of protein-ligand complexes can be better understood using RMSD plots, and it also sheds light on the stability and viability of these complexes as prospective therapeutic targets.Root mean square deviation RMSD plot of protein concerning ligands for 100ns MD simulation, showing most of the complexes are stable throughout the simulation trajectory, but compound 5533 showed less stability near 77-78 ns time frame; similarly, compound 4205 showed the same.The compound 4184 showed less stability at 55-56ns.The compounds 2771 and 4205, along with 5SE1, showed very stable at the end of the simulation as in the controlled 3152 compounds.The controlled compound 3152 showed a slight instability at 25-40 ns, but after 65 ns, it regained stability and maintained it till 100 ns.The compound 4205 shows a jump in stability at 75ns, and after that, it maintains that stability till 100ns.The RMSD plot for the rest of the protein-ligand complex showed stability, but the compound 5533 is very unstable throughout the trajectory and attained an RMSD value of 14-16 Å.The compound 4184-5SE1 showed stability throughout the trajectory with RMSD value ranging between 1-2 Å.All the compounds showed stability till 45 ns properly.After that, they attained conformational changes, which formed a stronger interaction between the protein and ligand and RMSD plot is shown as Fig 15.The root-mean-square-fluctuation (RMSF) is used to describe flexibility differences among residues.The backbone RMSF of each residue of the three systems, was calculated to analyse the backbone structure's flexibility.The vital amino acids between 10-25 show higher flexibility with an RMSF value of 0.5-0.6 nm.The average rmsf value is 0.35nm throughout the trajectory.The compound 2771 4205 showed the most negligible fluctuations throughout the trajectory compared to the controlled compound 3152.Amino acid residues of 5SE1: THR61, ILE62, Gly63, SER96, LEU97, GLN98, LYS99, LEU100, and PHE101 showed more flexibility in these two compounds 2160 and 4553 and their rmsf value ranging between 0.15-0.52nm.This indicated that the conformational changes of the systems at the end of simulations were relatively small.The RMSF is shown in Fig 16.A protein's volume and shape are described by the radius of gyration (Rg), which offers information on the compactness and stability of the protein's structure.A higher Rg value indicates an expanded system, which suggests a less dense and compact protein structure.A stable folded protein, on the other hand, usually shows a constant Rg value.Changes in the Rg value can signal structural changes in the protein, such as unfolding or denaturation, which can happen under various circumstances or due to various external factors.After analysing the plot, it was observed that the compactness of the compound 2771, 4205 is much less compared to the controlled 3152 with Rg value ranges between 3.2-3.75nm.This indicates that the compound 2771 and 4205 are better compact than 3152.The rest of the compounds showed The solvent-accessible surface area (SASA) plot analysis of the molecular dynamic simulation was used to gauge the protein's environmental exposure over time.The plot revealed information about the dynamics and structural changes of the protein during the simulation and is shown in Fig 18 .It was feasible to pinpoint areas of the protein with varying solvent exposure by keeping track of changes in SASA levels.This could point to protein-ligand interactions, structural rearrangements, or solvent-shielding effects.The SASA plot analysis helped to understand the behaviour of the protein better and offered vital data for researching its relationships and function in a dynamic context.In this study, after analysing solvent-accessible surface area for all the proteins and the compounds throughout the simulation, it was observed that all the systems have SASA values ranging between 53.96 nm2 and 302.11 nm2.Compounds 2771 and 4205 showed excellent solvent accessibility with 54-150 nm2 values.Overall, it can be concluded from this, which is comparatively reasonable compared to the controlled 3152 proteins, and the compound has excellent accessibility to the solvent environment.
The hydrogen bond interaction between the protein and ligand changes throughout the trajectory at 0ns, there are approx.1-2 h-bond, and the interaction increases near 17-18 ns due to some conformational changes of the compounds along with the conformational changes of the binding site residues.After 45ns, the interaction starts falling at 100ns; it maintained 2-3 interactions.The hydrogen bond interaction is shown in S1 Fig.The ligand interaction with the protein's binding site is shown in Figs 19-25.In all these images, it is observed that most of the compounds form h-bonds with the protein binding site and also show salt bridges.The compound 5533 also showed halogen bond interaction with the Ala135 residue of the protein's binding site.
Residue contact map and the mean smallest distance between C-alpha atoms of each amino acid residue 5SE1 corresponding to docked compounds.The residue contact map is analysed to calculate the smallest distance between c-alpha atoms that unify protein-ligand docked complexes that influence secondary structure elements to know their allosteric effects on the protein.The protein-ligand complex 5SE1-2771 5SE1-4205 showed atomic distance as shown in the above plots.From the above comparative analysis, it is evident that the interacting amino acids of the complexes were near form strong binding and stable conformations of the ligands towards the active pocket of 5SE1 and it is shown in Figs 26-32.The Free Energy Landscape (FEL) displayed an intricate energy landscape with several basins and transition states, pointing to a system with a rich conformational diversity.Numerous shallow wells for the compounds with protein signify metastable states, whereas substantial barriers signify transitions that require much energy.Notably, the conformation with the lowest energy is the most stable.According to this FEL study, the system investigates several structural configurations linked to a different energy level.Understanding the system's behaviour, including possible binding paths, conformational changes, and the thermodynamics driving molecular dynamics, is essential for comprehending biological processes and logical drug design.These insights into the energy landscape give vital information about the system's behaviour.After analysing the FEL plot, it was observed that the FEL achieved global minima of C-alpha backbone atoms of proteins concerning RMSD and rGyr.The FEL forecast for Principal Component investigation (PCA) was used to examine the MD simulations of the 5SE1 protein and certain chemicals' trajectories.The main aim was to understand atoms' complex, erratic motions in amino acid residues.This procedure revealed instances when the protein structure suffered deformation due to greater flexibility, shedding light on the variety in trajectories.This was accomplished by tracking the motion of internal coordinates in a threedimensional space, recording data throughout 100 nanoseconds, and presenting the data in a scatter pattern.Later, the idea of orthogonal sets or Eigenvectors was used to understand the controlled motion of each trajectory.In particular, the PCA analysis revealed an intriguing After analysing all the molecular docking results using PyMOL, MD Simulation trajectory analysis, it can be concluded that the citalopram and mirtazapine showed the best results in all aspects as we compared them with control 3152.It was observed that during simulation for 100ns, binding of these compounds to the active site of protein caused some conformational changes of the residues of the protein, due to which, from a reported time frame, it obtained stability till 100ns and binds appropriately as compared to earlier.

Molecular Mechanics Generalised Born Surface Area (MMGBSA) calculations
Further protein-ligand binding affinity was validated by calculating the complexes' MMGBSA binding free energies.MM/GBSA is emerging as a valuable and practical approach to predicting the binding energy.The Prime MM/GBSA module of the Schrodinger Suite was used to calculate the binding free energy of receptor-ligand complexes.After analysing the MMGBSA plot, it was observed that for most of the protein-ligand complexes, the values ranged between -40-50 region as shown in orange colour in the plot.Based on the MMGBSA result we can observe that the binding value for both 2771 and 4205 is between -30-40 which is very similar to the 5SE1-3152 controlled compound whose MMGBSA binding value also ranged between

Discussion
Molecular docking is a widely recognised technique in the field of structure-based drug design, enabling the examination of interactions between tiny molecules and molecular targets.By employing this methodology, a more comprehensive comprehension of the behaviour of small molecules within certain protein targets can be achieved, facilitating the ability to forecast the structure of ligand-receptor complexes.This capability holds significant importance in the realm of pharmaceutical research.The repurposing of pharmaceuticals has consistently demonstrated efficacy in the treatment of many medical issues and illnesses.SSRIs have demonstrated favourable outcomes in the treatment of AD.The present study employed molecular docking and subsequent molecular dynamics simulations to assess the interaction between SSRI medications and MARK4 protein.The objective was to identify the most favourable SSRI drug candidate capable of selectively targeting the overexpressed MARK4 protein in AD.The MMGBSA method was employed to calculate the binding energy of the top six medications, taking into account their binding affinity.This approach was utilised to verify the docking results.The results of the analysis indicated that these pharmaceuticals consistently and reliably interacted with the active site residues of MARK4 over the whole simulation trajectory.Mirtazapine, a pharmaceutical compound that specifically targets the MARK4 enzyme, exhibits potential for further exploration in various research settings, including in-vitro, in-vivo, and clinical investigations pertaining to AD.Previous studies conducted by Shamsi et al. in 2020 have demonstrated the binding affinity of donepezil, rivastigmine tartrate, and other cholinergic inhibitors to the active site of the MARK4 enzyme [27].According to a recent study conducted by Waseem et al. in 2021, it has been observed that the protein irisin, which plays a critical role in the decline of memory associated with AD, has a binding affinity towards MARK4 and contributes to the preservation of its stability.The user has provided a numerical reference [28].According to recent investigations [29], the molecules ganoderic acid A and ganoderenic acid B, produced from Ganoderma lucidium, have shown significant potential in their interaction with MARK4.According to the latest investigation, galantamine has been identified as a potentially effective inhibitor for MARK 4 [30].Based on our investigation, bolstered by previous research, it has been determined that mirtazapine, which exhibits affinity for MARK4, is regarded as the most auspicious SSRI medication among the 24 SSRI medications under consideration.Prior to conducting clinical trials, it is imperative to validate our findings through in vitro and in vivo research.

Conclusion
Currently, there are just five medications for AD that have been approved and SSRIs are used as off-label drugs in AD.In this study we presents a novel prediction that citalopram and mirtazapine are the most promising SSRIs from the list of 24 drugs which is targeting the MARK4 protein, which is related with AD.However, further validation is required through rigorous in vitro and in vivo investigations.